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OR 
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(3d or (three adj2 dimensional)) and 
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1 Multlresolution green's function methods for interactive simulation of large'Scale |j[ 
eiastostatic objects 
Doug L. James, Dinesh K. Pal 

January 2003 ACM Transactions on Graphics (TOG), volume 22 issue 1 

Full text available: ^ pdf(8.69 MB) Additional Information: full citation , abstract , references , citings , index terms 

We present a framework for low-latency interactive simulation of linear eiastostatic models, 
and other systems arising from linear elliptic partial differential equations, which makes it 
feasible to interactively simulate large-scale physical models. The deformation of the models 
is described using precomputed Green's functions (GFs), and runtime boundary value 
problems (BVPs) are solved using existing Capacitance Matrix Algorithms (CMAs). 
Multlresolution techniques are introduced to control the ... 

Keywords: Capacitance matrix, Green's function, deformation, eiastostatic, fast summation, 
force feedback, interactive real-time applications, lifting scheme, real-time, updating, 
wavelets 



2 Dynamics: Interactive physically based solid dynamics 
M. Hauth, J. GroB, W. StraBer 

July 2003 Proceedings of the 2003 ACM SXGGRAPH/Eurographics Symposium on 
Computer animation 

Full text available: ^ pdf(5.Q2 MB) Additional Information: full citation , abstract , references 

The interactive simulation of deformable solids has become a major working area in 
Computer Graphics. We present a sophisticated material law, better suited for dynamical 
computations than the standard approaches. As an important example, it is employed to 
reproduce measured material data from biologlcaf soft tissue. We embed it into a state-of- 
the-art finite element setting employing an adaptive basis. For time integration the use of an 
explicit stabilized Runge-Kutta method is proposed. 

3 Dynamic real-time deformations using space & time adaptive sampling 
Gilles Debunne, Mathieu Desbrun, Marie-Paule Cani, Alan H. Barr 

August 2001 Proceedings of the 28th annual conference on Computer graphics and 
interactive techniques 

Full text available: ^ pdf(1.45 MB) Additional Information: full citation , abstract , references , citin gs , index terms 

This paper presents a robust, adaptive method for animating dynamic visco-elastic 
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deformable objects that provides a guaranteed frame rate. Our approach uses a novel 
automatic space and time adaptive level of detail technique, in combination with a large- 
displacement (Green) strain tensor formulation. The body is partitioned in a non-nested 
multiresolution hierarchy of tetrahedral meshes. The local resolution Is determined by a 
quality condition that Indicates where and when the ... 

4 Collision s and def ormations: A multiresolution frannework for dynamic deformations 
Steve Capell, Seth Green, Brian Curless, Tom Duchamp, Zoran Popovic 
July 2002 Proceedings of the 2002 ACM SXGGRAPH/Eurographics symposium on 
Computer animation 

Full text available: ^ pdf(1.34 MB) Additional Information: full citation, abstract , references , citings, index temis 

We present a novel framework for the dynamic simulation of elastic deformable solids. Our 
approach combines classical finite element methodology with a multiresolution subdivision 
framework in order to produce fast, easy to use, and realistic animations. We represent 
deformations using a hierarchical basis constructed using volumetric subdivision. The 
subdivision framework provides topological flexibility and the hierarchical basis allows the 
simulation to add detail where It Is needed. Since vo ... 

Keywords: animation, deformation, hierarchical basis, multiresolution, physically-based 
animation, physically-based modeling 
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A study of level-of-detail in haptic rendering 
Shahram Payandeh, John Dill, Jian Zhang 

January 2005 ACM Transactions on Applied Perception (TAP), Volume 2 issue i 

Full text available: ^ pdf(1.03 MB) Additional Information: full citation , abstract , references , index tenms 

This paper presents an initial study of an approach to reduce computational overhead in 
haptic rendering of physically based models. Haptic rendering refers to the notion of adding 
physical properties and behavior, specifically a sense of touch or force feedback, to models 
of objects. In this way, a user through a haptic feedback device can feel interaction forces 
while visually observing the objects. Physically based modeling is particularly important 
when representing deformable objects. In th ... 



Keywords: Collision detection, deformable objects, haptic feedback, level-of-detail 



2 A scalable force pro pagation approach for web-based deformable simulation of soft 
tissues 

K. S. Choi, H. Sun, P. A. Heng, J. C. Y. Cheng 

February 2002 Proceeding of the seventh international conference on 3D Web 
technology 

Full text available: ^pdf(912.05 KB) Additional Information: full citation , abstract, references , index terms 

Physically based models and simulation are usually computationally intensive and not 
suitable for real-time interactive virtual reality applications including on-line medical 
training and surgical simulation. In this paper, we propose and develop a web-based 
scalable deformable model by simulating deformation of soft tissues as a successive force 
propagation process. This approach avoids laborious formulation of stiffness matrices in 
conventional mass-spring models. Computational speed is optimi ... 

Keywords: Java, Java 3D, deformable modeling, force propagation, scalability, surgical 
simulation 



3 Interaction and VR: A model for efficient and accurate interaction with elastic objects in Q 
ha ptic virtual environments 
Dan C. Popescu, Michael Compton 

February 2003 Proceedings of the 1st international conference on Computer graphics 
and interactive techniques in Australasia and South East Asia 
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Full text available: ^ pdf(260.21 KB) Additional Information: full citation , abstract, references 

This paper describes a method of modelling real-time interactions with elastic 3D objects 
represented by finite element models, which is particularly suitable for haptic virtual 
environments. The assumption we make is that the area of interaction of the external 
forces on the object is small. Our method provides a physically based solution and only 
requires the precomputation of the inverse of the stiffness matrix. It can be naturally 
coupled with a technique of local multiresolution collision d ... 

Keywords: deformable objects, finite element method, haptics, linear elasticity, virtual 



4 Multiresolution green's function methods for interactive simulation of large-scale 
elastostatic ob j ects 
Doug L. James, Dinesh K. Pal 

January 2003 ACM Transactions on Graphics (TOG), volume 22 issue 1 



We present a framework for low-latency interactive simulation of linear elastostatic models, 
and other systems arising from linear elliptic partial differential equations, which makes it 
feasible to interactively simulate large-scale physical models. The deformation of the 
models is described using precomputed Green's functions (GFs), and runtime boundary 
value problems (EVPs) are solved using existing Capacitance Matrix Algorithms (CMAs). 
Multiresolution techniques are introduced to control the ... 

Keywords: Capacitance matrix, Green's function, deformation, elastostatic, fast 
summation, force feedback, interactive real-time applications, lifting scheme, real-time, 
updating, wavelets 



^ Applicati ons: Collision detection and tissue modelin g in a VR-simulator for eye sur g er y Q 
Clemens Wagner, Markus A. Schill, Reinhard Manner 

May 2002 Proceedings of the workshop on Virtual environments 2002 

Full text available: ^ pdf(3>43 M B ) Additional Information: full citation , abstract , reference s, index terms 

This paper gives a survey of techniques for tissue interaction and discusses their application 
In the context of the Intra-ocular training system EyeSI. As key interaction techniques 
collision detection and soft tissue modeling are identified. For collision detection In EyeSi, an 
enhanced image-based approach for collisions between deformable surfaces and rigid 
objects is presented. By exploiting the computing power of graphics processing units, it 
achieves higher performance than existing geome ... 

^ 1-2 VRC in bio & medical sciences: A haptic needle manipulation simulator for Chinese Q 
acupuncture learning and training 

Pheng-Ann Heng, Tien-Tsin Wong, Ka-Man Leung, Yim-Pan Chui, Hanqlu Sun 
June 2004 Proceedings of the 2004 ACM SIGGRAPH international conference on Virtual 
Reality continuum and its applications In industry 

Full text available: " ^DdfOSe-IS KB) Additional Information: full citation , abstract , references , index terms 

This paper presents a haptic needle manipulation simulator for Chinese acupuncture 
learning and training. Students can learn and practise acupuncture in the proposed 3D 
interactive virtual environment that supports force feedback interface for needle insertion. 
So that, students not only "see" but also "touch" the virtual patient. With the high 
performance computers, highly informative and flexible visualization of acupuncture points 
of various related meridian and collateral can be highlighted ... 
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7 Modeling and animating hands & bodies: Construction and animation of anatomically 
based human hand models 
Irene Albrecht, Jorg Haber, Hans-Peter Seidel 

July 2003 Proceedings of the 2003 ACM SIGGRAPH/Eurographics Symposium on 
Computer animation 

Full text available* 153 Ddf(7 55 MB) Additional Information: full citation , abstract , references , citings , index 
' ^ terms 

The human hand Is a masterpiece of mechanical complexity, able to perform fine motor 
manipulations and powerful work alike. Designing an animatable human hand model that 
features the abilities of the archetype created by Nature requires a great deal of anatomical 
detail to be modeled. In this paper, we present a human hand model with underlying 
anatomical structure. Animation of the hand model is controlled by muscle contraction 
values. We employ a physically based hybrid muscle model to convert ... 

^ Ha ptic rendering: programmin g touch interaction with virtual objects 
K. Salisbury, D. Brock, T. Massie, N. Swarup, C. Zilles 

April 1995 Proceedings of the 1995 symposium on Interactive 3D grapliics 

Full text available: Ifi p df (897.60 KB) Additional Information: full citation, abstract, references , citings, index 

terms 

Haptic rendering is the process of computing and generating forces in response to user 
interactions with virtual objects. Recent efforts by our team at MIT's AI laboratory have 
resulted in the development of haptic interface devices and algorithms for generating the 
forces of interaction with virtual objects. This paper focuses on the software techniques 
needed to generate sensations of contact interaction and material properties. In particular, 
the techniques we describe are appropriate fo ... 

^ Dynamics: Interactive physically based solid dynamics 
M. Hauth, J. GroB, W. StraBer 

July 2003 Proceedings of the 2003 ACI^ SXGGRAPH/Eurograpliics Symposium on 
Computer animation 

Full text available: ^ pdf(5.02 MB) Additional Information: full citation , abstract , references 

The Interactive simulation of deformable solids has become a major working area in 
Computer Graphics. We present a sophisticated material law, better suited for dynamical 
computations than the standard approaches. As an important example, It Is employed to 
reproduce measured material data from biological soft tissue. We embed it Into a state-of- 
the-art finite element setting employing an adaptive basis. For time integration the use of 
an explicit stabilized Runge-Kutta method is proposed. 

Dynamic real-time deformations using space & time adaptive sampling 
Gilles Debunne, Mathieu Desbrun, Marie-Paule Cani, Alan H. Barr 

August 2001 Proceedings of the 28th annual conference on Computer graphics and 
interactive techniques 

Full text available: fi ^pdf(1.45 MB) Additional Infomriation: full citation, abstract , references , citings, index 
^ ^^^^^ 

This paper presents a robust, adaptive method for animating dynamic visco-elastic 
deformable objects that provides a guaranteed frame rate. Our approach uses a novel 
automatic space and time adaptive level of detail technique, in combination with a large- 
displacement (Green) strain tensor formulation. The body is partitioned in a non-nested 
multiresolution hierarchy of tetrahedral meshes. The local resolution is determined by a 
quality condition that indicates where and when the ... 

''I Modellin g : Bag-of-particles as a deformable model 
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D. Stahl, N. Ezquerra, G. Turk 

May 2002 Proceedings of the symposium on Data Visualisation 2002 

Full text available: ^ Ddff6.32MB^ Additional Information: Ration, abstriet. references , citings, indix 

We present an interactive, physically-based, elastlcally deformable model using a particle 
systenn to model surfaces with interior volumes that can be haptlcally felt. Oriented 
particles used in existing surface-only models, and unoriented particles used in volume-only 
simulations are combined to form a bag-of-partlcles. Multiple species of surface and volume 
particles, coupled with prede£ned interspecies parameters, determine the elastic properties 
of a bag. Starting with an object represe ... 

^2 Cliaracter/web: Animated deformations with radial basis functions 

Jun-yong Noh, Douglas Fidaleo, Ulrich Neumann 

October 2000 Proceedings of the ACM symposium on Virtual reality software and 
technology 

Full text available: ^ Ddf(2.26 MB) Additional Information: full citation , abstract , references , index terms 

We present a novel approach to creating deformations of polygonal models using Radial 
Basis Functions (RBFs) to produce localized real-time deformations. Radial Basis Functions 
assume surface smoothness as a minimal constraint and animations produce smooth 
displacements of affected vertices in a model. Animations are produced by controlling an 
arbitrary sparse set of control points defined on or near the surface of the model. The 
ability to directly manipulate a facial surface with a small numbe ... 

Keywords: Facial Animation, Geometry Deformation, MPEG-4, Radial Basis Functions 



13 Simulation and preferences: Real-time simulation of elastic objects in virtual 
environments using finite element method and precomputed Green's functions 
Igor NIkitin, Lialia Nikitina, Pavel Frolov, Gemot Goebbels, Martin Gobel, Stanislav Klimenko, 
Gregory M. Nielson 

May 2002 Proceedings of the workshop on Virtual environments 2002 

Full text available: ^ pdf(1.26 MB) Additional Infonnation: full citation , abstract , references , citings 

Simulation of an object's elastic deformation is an important feature in applications where 
three-dimensional object behavior is explored. In addition, the benefits of user-object 
interactions are best realized in interactive environments which require the rapid 
computation of deformations. In this paper we present a prototype of a system for the 
simulation of elastic objects in Virtual Environments (VE) under real-time conditions. The 
approach makes use of the method of finite elements and prec ... 
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We present a novel framework for the dynamic simulation of elastic deformable solids. Our 
approach combines classical finite element methodology with a multiresolution subdivision 
framework in order to produce fast, easy to use, and realistic animations. We represent 
deformations using a hierarchical basis constructed using volumetric subdivision. The 
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subdivision framework provides topological flexibility and the hierarchical basis allows the 
simulation to add detail where it Is needed. Since vo ... 
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Clothing is a fundannental part of a character's persona, a key storytelling tool used to 
convey an intended impression to the audience. Draping, folding, wrinkling, stretching, etc. 
all convey meaning, and thus each is carefully controlled when filming live actors. When 
making films with computer simulated cloth, these subtle but important elements must be 
captured. In this paper we present several methods essential to matching the behavior and 
look of clothing worn by digital stand-ins to their ... 
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June 2004 Proceedings of the 2004 ACM SIGGRAPH international conference on Virtual 
Reality continuum and its applications in industry 
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Realism of the mass-spring model based simulations remains questionable. In this paper, 
we develop the relationship between the material properties of a continuum plate and the 
values of the static parameters for its mass-spring model. The result indicates the necessity 
to use nonzero preload for the springs. Then we apply the linear approximation of the result 
to assign parameters for the mass-spring model of a rectangular plate. The simulation 
shows that the mass-spring system with this set of ... 
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We investigate techniques for modeling contact between quasi-rigid objects - solids that 
undergo modest deformation in the vicinity of a contact, while the overall object still 
preserves its basic shape. The quasi-rigid model combines the benefits of rigid body models 
for dynamic simulation and the benefits of deformable models for resolving contacts and 
producing visible deformations. We argue that point cloud surface representations are 
advantageous for modeling rapidly varying, wide area c ... 
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This demonstration presents the current state of an on-going team project at Simon Fraser 
University in developing a virtual environment for helping to train surgeons in performing 
laparoscopic surgery. In collaboration with surgeons, an initial set of training procedures 
has been developed. Our goal has been to develop procedures in each of several general 
categories, such as basic hand-eye coordination, single-handed and bl-manual approaches 
and dexterous manipulation. The environment is based ... 

Keywords: haptics, surgery training, surgical simulation, virtual laparoscopy, virtual reality 
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Tissue Elasticity Reconstruction Based on 
Ultrasonic Displacement and Strain Images 

A. R, Skovoroda, S. Y. Emelianov, and M. O'Donnell, Fellow, IEEE 



Abstract — A method is presented to reconstruct the elastic 
modulus of soft (issue based on ultrasonic displacement and strain 
images. Incompressible and compressible media are considered 
separately. Problems arising with this method, as well as ap- 
plications to real measurements on gel-based, tissue equivalent 
phantoms, are given. Results show that artifacts present in strain 
images can be greatly reduced using a hybrid reconstruction 
procedure based on numerical solution of the partial dififeren- 
tial equations describing mechanical equilibrium of a deformed 
medium. 



I. Introduction 

IMAGING internal soft tissue displacements and strains re- 
sulting from mechanical forces applied to the body surface 
is rapidly developing into a new diagnostic modality [1H32]. 
Internal deformational images, however, emphasize both the 
spatial distribution of the Young*s or shear modulus and 
global boundary conditions, including mechanical constraint 
of the body, its geometry, the types of external and internal 
forces, etc. That is, displacement and/or strain images may 
exhibit significant artifacts due to global boundary conditions, 
as discussed in [28], [29]. The primary objective of the 
work described here is to significantly reduce artifacts in 
elasticity images by directly reconstructing and imaging the 
elastic Young's modulus. Although absolute quantitation may 
ultimately be important for certain applications, differential 
diagnosis based on tissue elasticity will probably be based 
on relative modulus reconstruction. We strongly believe that 
artifacts due solely to global boundary conditions must be min- 
imized for elasticity imaging to become a routine diagnostic 
procedure. 

Computing the mechanical properties of a medium based 
on its response to mechanical action can be posed in a 
number of ways. In this paper we have used two different 
formulations of elasticity reconstruction. The first starts with 
the complete spatial distribution of the strain tensor or dis- 
placement; the second uses limited experimental data such, as 
an image of a single strain component. Because of physical 
limitations inherent in measuring internal displacements and 
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strains with an ultrasound imaging system, i.e. limitations 
of traditional longitudinal speckle tracking algorithms for 
large absolute displacements and the poor accuracy of lateral 
displacement measurements (due to lower spatial frequencies 
laterally), worie to date has focused on estimating elastic 
moduli with limited deformational infonnation [17], [21], 
[24], [28], [29]. Based on simpUfied models of both the 
elastic modulus distribution in the body and the mechanical 
boundary conditions, diese methods are accurate ozily in 
limited applications or, otherwise, produce laige artifacts in 
the elasticity reconstruction [29]. Consequently, expanded 
reconstruction methods are needed to handle more complex 
objects and boundary conditions. In this paper, we explore 
new reconstruction methods based on quantitative strain and 
displacement images. 

In the dieory section, the general approach to reconstruction 
based on the conunon model of linear, elastic, isotropic, 
incompressible media is presented Compressible media are 
considered in Appendix I. Following that, a hybrid reconstruct 
tion procedure independent of global boundary conditions is 
described. Also presented here are practical methods to esti- 
mate the elastic modulus using only limited experimental data, 
i.e., only accurate measurements of the axial displacement 
and longitudinal strain. Some additional cases of elasticity 
reconstruction based on limited measurements are considered 
in Appendix n. The details of the hybrid reconstruction method 
based on internal tissue displacement and strains estimated 
from ultrasonic images is described in the next section. In the 
results section, a specific reconstruction algorithm appropriate 
for a plane strain state is tested using experimental images of 
internal deformations in tissue equivalent, gel-based phantoms. 
These displacement and strain images were made with an 
ultrasound-based deformational imaging system described in 
[27]-[29]. The paper concludes with a discussion of the results. 

IT. Theory 

Consider a three-dimensional (3-D) volume V of deformed 
media with the displacement vector U = (ui, U2, ws) in 
Cartesian coordinates X = (zi, X2, xa). Volume V can be 
eidier the entire mechanical body, or a region of interest (ROI) 
inside the object under study. 

The most general form of Newton's 2nd law describing the 
motion of a mechanical body under static deformation (i.e., 
the equilibrium condition) is, 

3 

X)^0'>+A^O. 1 = 1,2,3 (1) 
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where aij is one component of the 2nd ranked stress tensor and 
fi is the body force per unit volume acting on the body in the 
Xi direction [29], [33]-[351. This model applies for soft tissue 
at a spatial scale sampled by diagnostic ultrasound (i.e.. at a 
scale greater than or comparable to an ultrasound wavelength). 
In (1), and the entire paper, the lower index after a conmia 
means differentiation with respect to the coiresponding spatial 
coordinate. Equation (1) must be satisfied at every internal 
point of the body. We should also note here that most soft 
tissues and tissue-like materials (with the exception of lung, 
for example) can be considered incompressible [36], i.e., the 
Poi$son*s ration is equal to 0.5. In the main body of the 
paper only incompressible media are considered. Compress- 
ible media (i.e., Poisson's ratio less dian 0.5) are discussed 
separately in Appendix I. Treating slightly compressible media 
as completely incompressible greatly simplifies the partial 
differential equations describing internal deformations and 
permits stable numerical solution of these equations. A com- 
plete compressible description, although theoretically more 
accurate, often leads to unstable numerical solution in the limit 
that the Poisson's ratio approaches 0.5. 

Assuming linear elasticity, the components of the stress 
tensor in an isotropic, continuous, incompressible medium 
under static deformation are: 

<^ii =P^ij + 2^ij, (2) 

where p is the static internal pressure, defined as [37] 

Jini^lA(dit;U)] = p. (3) 

dttiu-*o 

In (2) 6ij is the Kronecker delu symbol and Cij is one 
component of the 2nd ranked symmetric strain tensor, defined 
as 

where u, is one component of displacement vector U = 
(ui, U2, tx^). Note that a pressure term must be included in (2) 
to fully describe deformations in incompressible media. Also, 
some tissues, such as muscle, are anisotropic [35], necessitat- 
ing a more general stress-strain relationship [33], [34]. In this 
paper, however, only isotropic media are considered. 

In (2) and (3) the parameters A and are Lame coefficients. 
In general, the longitudinal Lame coefficient A and shear 
modulus /i describe the elastic behavior of a mechanical 
body. However, with assumption (3), any statically deformed, 
isotropic, continuous, incompressible mechanical body can 
be completely characterized by a spatial distribution of a 
sin^e material parameter, either shear modulus /i or Young* s 
modulus E since they are simply proportional to each other, 
i.e. £ = 3^. Note that this is a fundamental difference 
between incompressible and compressible media. That is, 
the Young's modulus completely describes the static elastic 
properties of soft tissue, where its value may vary widely 
between different tissues (see, for example, reference [31]). 



In addition, the incompressibility condition leads to the 
following relationship between displacement or strain com- 
ponents 

divU = £n + £22 + ^33 = +U2,2 +«3i3 = 0. (5) 

The main goal of elasticity imaging is to reconstruct the 
elastic modulus of any desired tissue region using available 
measurements of strain and displacement components. The 
Young*s modulus £ is an arbitrary function of position, i.e. 

£;(X) = Eok(X) (6) 

where £^ is a constant and A:(X) is not generally a continuous 
function. Elasticity image reconstruction, therefore, centers on 
estimating k{X). 

Elasticity reconstruction in an isotropic medium must ac- 
curately represent both continuous and discontinuous changes 
in the modulus. Both types of spatial variation are illustrated 
in Fig. 1 for a single abnormality (i.e., inclusion) in an 
otherwise homogeneous medium. Small fluctuations in the 
modulus throughout the image plane represent local changes 
not indicative of true material differences. In Fig. 1(a), a 
clearly bounded inhomogeneity ("inclusion") with a discon- 
tinuous change in elastic modulus is illustrated. For this 
type of abnormality the magnitude of spatial changes in the 
elastic modulus at the boundary is much larger than any 
small scale fluctuations present. In contrast, a continuous 
inclusion with finite derivatives everywhere is presented in 
Fig. 1(b). Depending on the size and elasticity contrast of 
the inclusion, derivatives in the elastic modulus may be only 
slightly larger than small scale fluctuations in the medium. 
Both types of inclusion are important, however, and represent 
reasonable assumptions for a large set of tissue abnormalities. 
Due to elasticity variability and errors in strain measurement, 
reconstruction approaches for these two types of abnormalities 
may differ. 

A. Bounded Inhomogeneities 

We focus first on local, deariy bounded inclusions having 
a boundary G and residing in tissue with otherwise smoothly 
varying mechanical characteristics, as illustrated in Fig. 1(a). 
For this type of inclusion the unknown A:(X) is not a continu- 
ous function, i.e. it has a discontinuity at the boundary G. Note 
that the position of the boundary G is also unknown a priori 
and needs to be determined. Therefore, the reconstruction 
problem reduces to detecting the unknown inclusion bound- 
ary and, consequently, reconstructing the unknown Young^s 
modulus at this boundary. 

The stress continuity condition at the boundary G has the 
form [34]: 

Z WiMi = E K- * - = 0, < = 1, 2, 3 (7) 

where the square parentheses denote discontinuity of the 
corresponding external and internal terms at boundary points 



SKOVORODA a at^ TISSUE ELASTICITY RECONSTRUCTION BASED ON ULTRASONIC DISPLACEMENT AND SITIAIN IMAGES 



749 




Xg C G, and rij is the jth component of the unit normal 
vector n = (ni , n2 , na) at the boundary. In this and the 
following equations superscripts "mt*' and "exf * refer to 
internal and external variables at boundary points. Note that for 
any other point of the volume V inside or outside the inclusion, 
the stress continuity condition (7) is automatically satisfied for 
any arbitrarily chosen direction of a normal vector. 

Assuming an incompressible, isotropic medium, substituting 
expression (2) into (7), and combining the results, the stress 
continuity condition becomes: 

= P^'^na + 2 /i'"'(eS*ni + eg^nj + e^'m). (8) 

Since the pressure p cannot be directly measured by imaging 
systems, it is analytically eliminated from the stress disconti- 
nuity equations resulting in the following pair of equations 

[26]: 

r[nin2(e;f -e|f) + [(n2)2-(n02] 

• eff + n3(n2e5f - metS')! = nmaCeir - 4?) 
+ [(^2)' - (ni)^]er2* + n3(n2£l^' - meg') 

r[nin3(ejf ^£5f) + [(n3)2-(ni)2j 

• ^11' + nzCnaetr - nie|f )] = ninz{e[^' - 4?) 

[(na)^ - (ni)2]ei5* + njCna^?* - niej?). (9) 

In (9) the ratio of shear or Young*s moduli, = 
£€xt^£%nt^ is denoted by F. Compressible media are treated 
separately in Appendix I. 

Both inclusion boundary detection and Young's modulus 
reconstruction at this boundary are based on (9) and performed 
in two steps. First, experimentally derived (9) are solved with 
respect to F for an arbitrarily chosen direction of the normal 
vector n = (ni, ns, TI3). At any boundary points where the 
direction of the **testing** normal vector coincides with the 
direction of the real normal vector, the true value of T is 
obtained. At other points, however, the value F will not rep- 
resent the exact value but will highlight mclusion boundaries. 



This procedure can be repeated for several directions of the 
normal vector to complete boundary detection. Next, after 
the inclusion boundary is identified, the true value of P can 
be calculated at boundary points since the direction of the 
normal vector will be predicted after the boundary detection 
procedure previously described. Note that if all components 
of the strain tensor are known, the first equation in (9) can 
be used to detect the boundary of the inclusion and then 
determine the value of F (i.e., relative Young's modulus ratio 
at this boundary), and the second one to estimate the accuracy 
of experimental strain measurements in the neighborhood of 
the boundary point Xg- Also, if the mechanical properties 
do not vary significantly inside the inclusion, reconstruction is 
complete since the Young's modulus at the boundary is found. 

Equation (9) can also be used for boundary detection of 
more complicated elasticity inhomogeneities. If the Young's 
modulus changes smoothly both inside and outside the inclu- 
sion, then experimentally derived (9) solved with respect to F 
for any chosen direction of the vector n will give smoothly 
varying results everywhere except at boundary intersection 
points. Although the value F will not represent the exact 
value, this will highlight inclusion boundaries. In general, this 
approach produces a qualitative reconstruction, i.e., boundary 
detection, and will require Young's modulus reconstruction at 
every point of the volume V or any chosen region of interest 
This is considered in the following section. Nevertheless, 
such information derived from measured strain images can be 
extremely inqx>rtant in those cases where the inclusion or other 
tissue structures cannot be directly visualized by the imaging 
system [28], [29]. 

In (8), (9) a linear stress-strain relation is assumed, but 
no assumption is made about the strain-displacement relation. 
Also, only the strain distribution in the immediate neighbor- 
hood of the boundary is needed to reconstruct the Young's 
modulus at this boundary. Consequently, this method can be 
used without requiring small strain magnitude as long as the 
stress-strain relation is still approximately linear, i.e., described 
by (2). In most soft elastic materials such as tissue, the strain 
reaches a limiting value ed> where the strain-displacement 
relation (4) is no longer valid, before the stress-strain relation 
goes nonlinear. Only upon reaching a second, higher limiting 
strain value, £«, do both single strain-displacement (4) and 
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linear stress-strain (2) relations breakdown [33], [35], [38], 
[39]. Therefore, methods based on stress continuity (8), (9) 
are valid up to relatively high strain magnitude (for example, 
up to 20% [38]). Note that limiting strain values vary for 
different tissues. The use of high strains is extremely important 
in clinical applications where laige displacements can enhance 
the signal to noise ratio of quantitative strain images [28]. 

Equation (9) shows that for an arbitrary deformed state, 
reconstructing boundary values of the Young*s modulus for 
clearly bounded inhomogeneities requires measurement of all 
strain components in the neighborhood of the boundary G. For 
a general, three-dimensional inclusion, this means complete 
three-dimensional strain data. To reconstruct the Young's 
modulus in a single plane, all displacement components must 
only be measured in three close parallel planes intersecting 
the volume V. From measurements in these three planes, 
coefficients needed for reconstruction in the central plane 
based on (9) can be calculated numerically with second 
order accuracy. General three-dimensional reconstruction of 
the inclusion boundary can be done by moving the imaging 
planes in parallel. Note tiiat only a pair of close parallel planes 
can also be used for numerical reconstruction in a single plane, 
but will, however, exhibit greater error since less accurate 
finite differences will be used to ^proximate all derivatives. 

B. Distributed Inhomogeneities 

Not all inhomogeneities will exhibit a clear boundary and, 
therefore, reconstruction of an isotropic mediimi with spatially 
continuous changes in the Young's modulus is required. An 
example of such a distribution is given in Fig. 1(b). For this 
type of reconstruction the function k{X) in (6) is assumed to 
be continuous with well defined spatial derivatives. 

Using the equations defining a linear, incompressible elastic 
medium, the equilibrium condition (1) can be rewritten in the 
form [26]: 

2ei2(fe,ll -fc,22 ) + 2(tl2,2 -1*1,1 )fe,12 

+ 2e23A;,i3 -2£i3fc,23 +(AU2 + Wi2,l )A:,i 

- (Atli - a;i2,2 )k,2 +Wi2,3 *;,3 -f Aa;x2fe + F12 = 0 

2€l3(fc»ll -^)33 ) + 2e23fe>12 

+ 2(u3,3 -uui )A;,i3 -2ffi2A:,23 

+ (Ali3 + Wi3,i)A;,i+Wl3,2A;,2 

- (Aui - wi3,3 )fe,3 +Au;i3A; + F13 = 0 

2€23(^i22 -k,33 ) + 2€iik,i2 -2ei2fc,13 

+ 2(W3,3 -1X2,2 )^,23 +'*'23,1 +(AU3 + LJ23,2 )fc,2 
~(A«2 "~ '*'23>3 )^}3 +Au;23^ + ^23 = 0 (10) 

where the pressure p is eliminated from this system of equa- 
tions, and the following notations are used: 

W\2 =U2,1 -Ui,2 ,ti^l3 = W3,l -"1,3 it»^23 = "3,2 -«2»3 » 

A =aV(ax?) + d^idxl) + 8^(8x1), 
Fa2 =(3/Eo){/2,i -/i,2),Fi3 = (3/£7o)(/3,i-/i,3), 

^23 = (3/^^o)(/3,2 -/2,3 )• 

Analytical solution of (10) is not generally possible for an 
arbitrary spatial distribution of the elastic modulus k{X), Note 



also that this set of equations represents a boundary value 
problem and has an infinite number of solutions unless some 
boundary conditions for k(X) are specified. 

Complete numerical solution of these equations is possible 
if all components of the displacement vector and strain tensor 
are measured. It also requires that the boundary conditions for 
k{X) be specified in the volume V. However, as considered in 
Appendix 11, for some cases of Young's modulus distribution, 
numerical solution of (10) can be obtained in a single plane 
by discretizing all functions and computing spatial derivatives 
with finite differences to yield a set of simultaneous linear 
equations. That is, the complete displacement vector U must 
be accurately measured in at least three close parallel planes 
intersecting the volume V. Given complete displacement data 
in these planes and boundary conditions for k{X) only in the 
central plane, all coefficients of the differential equations can 
be calculated numerically with second order accuracy and the 
elasticity distribution can be reconstructed in the central plane. 

C Hybrid Reconstruction 

In practice, neither clearly bounded inclusions nor smoodily 
distributed inhomogeneities are complete models of elastic 
modulus variations. To handle more general cases, a hybrid 
procedure is used. 

First, strain images are processed based on (9) to highlight 
boundaries between regions of different elastic modulus. If 
clearly bounded inclusions are present, then the Young's 
modulus can be obtained directiy from (9) at points along the 
boundary. More importantly, this procedure can define regions 
of very small modulus variations (i.e., where F is close to one). 

Following boundary detection, closed contours of small 
elasticity variations are defined. The total dimension of the 
contour must include several resolution cells of the final 
reconstruction. The modulus along these contours can be 
considered constant thereby providing a boundary condition 
for complete reconstruction of the elasticity within this region 
of interest (ROI) based on numerical solution of (10). The 
elasticity reconstructed in this way is the modulus relative 
to the modulus along the boundary. If the elastic modulus 
along this boundary is known (for example, if a specific tissue 
type such as fat constitutes the boundary region), then the 
reconstructed modulus is absolute. If the modulus along the 
boundary is not known, then the reconstruction remains rela- 
tive. Even if the reconstruction is relative, potential artifacts 
in pure strain images due to tiie details of global boundary 
conditions (i.e.. geometry of the object, types of external load, 
etc.) are removed. 

Solution of (10) requires that both first and second order 
spatial derivatives of A:(X) be finite. If the ROI defined using 
(9) contains a bounded inclusion, then this may not be satisfied. 
To avoid this, a continuous elasticity distribution is forced by 
low pass filtering both displacement and strain data within the 
ROI prior to numerically solving (10). This filter effectively 
smoothes the reconstruction so that k{X.) is forced to be 
continuous with finite first and second order spatial derivatives. 
A negative consequence is that the reconstructed modulus 
is underestimated both along the boundary of large, cleariy 
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bounded inclusions and within inclusions of spatial extent 
comparable to the dimensions of the filter. If necessary, the 
modulus can be estimated at these points from solution of (9), 
i.e., firom the boundary detection process. This hybrid proce- 
dure allows local reconstruction of the elastic modulus within 
the body without knowledge of global boundary conditions. 

D. Plane Strain State 

If all components of the displacement vector and strain 
tensor are available, then the Young's modulus can be recon- 
structed with the hybrid procedure outlined above. However, 
current strain imaging systems based on ultrasonic speckle 
tracking are two-dimensional [17], [28], [29], and conse- 
quently, boundary detection and Yoimg's modulus recon- 
struction within selected ROFs must be estimated based on 
displacement and strain data from a single imaging plane. 
To reconstruct the elastic modulus using current ultrasonic 
equipment either some symmetries must be assumed in the 
distribution of the Young's modulus, or the deformation pat- 
tern must be carefully controlled. Some simplifications to the 
general equations based on spatial symmetries in k{X) per- 
mitting reasonable reconstruction from only two-dimensional 
displacement and strain measurements are considered in Ap- 
pendix n. Below, we consider a specific type of deformation 
requiring displacement and strain data only in the (xi, 0) 
plane (i.e., = 0) to reconstruct the elastic modulus in that 
plane. 

In some practically important cases the out of plane dis- 
placement U3 is either zero or small compared to the others 
and the two in-plane components ui and U2 do not vary 
significantly as a function of the out-of-plane coordinate. This 
is possible, for example, if one dimension of the region far 
exceeds the odiers and the character of the inhomogeneity as 
well as the manner of loading does not change or changes 
only slightly along this direction. An example of such a 
deformation is given in Fig. 2. Here a uniform surface dis- 
placement is applied along the long xz axis of a cyhndrical 
medium. Fig. 2(a) shows a three-dimensional view of this 
deformation geometry and Fig. 2(b) and (c) illustrate two 
oithogonal cross sections. IVo different cases are presented: 
a cylindrical inhomogeneity extending the entire length of the 
object and a spherical inhomogeneity placed at the center. For 
the cylindrical inhomogeneity, this deformation will accurately 
approximate a plane strain state for imaging at the central plane 
X3 = 0 [28], [29]. For the spherical inclusion, it produces 
a weak approximation to a plane strain state. As will be 
demonstrated below, if a uniform surface displacement is 
applied then a plane strain state is a good approxiniation even 
if the inhomogeneity is not extremely long. 

For a plane strain state, in-plane components ui and U2 of 
the displacement vector U are functions only of x\ and X2, and 
U3 = 0. Also, at the inclusion boundary only two components 
of the normal vector are nonzero, i.e., n = (ni, n2, 0). For 
these conditions, the second equation of (9) is automatically 
satisfied and the first one leads to the following. 



In this case, (11) should be satisfied along the entire boundary 
of a clearly bounded inclusion contained in the plane. 

Some additional simplifications based on the condition of 
incompressibility (5) are possible. For a plane strain state, 
the condition of incompressibility reduces to en -h £22 = 0, 
and consequently, at any point along the boundary where two 
components of the normal vector are equal in magnitude, i.e., 
(ni)' = (n2)^. (11) reduces to 

^ elf ^^^^ 

The parameter F is defined above as the ratio of Young* s 
moduli of internal and external media at the boundary of the 
inclusion. That is. at any point where two components of the 
unit normal to the border of a clearly bounded inclusion are 
equal in magnitude, the Young* s modulus can be accurately 
computed from a single component of the strain tensor. This 
is a key result since only a single component is accurately 
measured with current strain imaging systems [17], [28]. 

Accurate reconstruction of the Young's modulus based on 
(12) is only possible at special points. However. (12) can be 
generally used to detect inclusions since there will be large 
changes in V near the boundary over a wide range of n\ 
and 712. Following boundary detection, all points satisfying 
(ni)2 = {n^f can be identified. Al these points the Young's 
modulus obtained from (12) represents an accurate estimate at 
the boundary. For all other points, the shear strain component 
€12 should be taken into accounL Nevertheless, the simple 
expression presented in (12) can be used to detect inclusion 
boundaries even if the inclusion exhibits significant but smooth 
spatial changes in the Young's modulus. 

Following boundary detection and definition of a closed 
curve of constant elastic modulus, (10) must be solved nu- 
merically within all ROFs. For a plane strain state, where 
the displacement vector components are = ui(a:i, X2), 
U2 = U2{xu 3C2), U3 = 0, the set of (10) reduces to a single 
nontrivial equation: 

{kai - k,22 )(mi,2 +W2,i ) + 2A;,i2 (u2,2 ) 

+ 2k,i (ti2,ll +tA2,22 ) - 2fc,2 (Ui,ii -HUi,22 ) 

+ k{u2an -f «2,22i -«i,U2 -^1,222 ) 
+ Fi2 = 0 

P23=0. (13) 

Note that only two-dimensional in-plane displacement and 
strain measurements are needed to reconstruct the modulus in 
the plane 0:3 = 0. 

in actual ulti^onic strain imaging systems, one component 
of the displacement is measured with much higher accuracy 
than the other. This can create reconstruction error based 
on (13). Fortunately, this problem can be largely overcome 
because the incompressibility condition for a plane strain state 
is simply ui,i + U2.2 = 0, where U2,2 is the longitudinal strain 
component accurately measured within the image plane. Based 
on this, accurate measurements of the lateral ui component of 
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<C) 

Fig. 2. Schematic of a uniform dcfonnation applied along the long X3 axis of a cylinder, (a) Three-dimensional view of this deformation geometry. Two 
different cases are presented: a cylindrical inbomogencity extending the entire length of the object, and a spherical inhomogencity placed at the center, where 
a plane strain state is well approximated at the central cross-section of the object (b) and (c) Illustrate two orthogonal cross-sections through the same object. 



the displacement vector are needed along only one longitudi- 
nal curve to compute the in-plane distribution of the lateral 
displacement £rom the incompiessibility equation. 

The specific approach used to solve the 2nd order partial 
differential in (13) depends on the relative magnitudes of 
the coefficients in [40]. It is generally hypert)olic, although 
a degenerate (i.e., parabolic) case is possible. For an incom- 
pressible material, the degenerate case for the plane strain 
state corresponds to an in-plane translation with rotation of 
the volume as a rigid body. This means that the distribution 
of relative Young's modulus k{X) cannot be reconstructed 
without knowledge of the internal pressure distribution for 
this special case, and, therefore, either the pressure must be 
measured and included into reconstruction, or different strain 
and displacement patterns must be created within the object 
using different boundary conditions. Note thiat there are some 
very specific cases where (13) will always be degenerate 
regardless of boundary conditions and, therefore, will require 



measurements of the internal pressure distribution. Generally, 
however, the boundary defined from (12) produces a nonde- 
generate boundary value problem with a unique solution. Note 
also that (13) for the nondegenerate case r^resent a boundary 
value problem and, therefore, boundary conditions for ^(X) 
must be specified to determine the unique solution. 

III. Methods 

Experiments were performed on a number of gel-based 
cylindrical phantoms; a detailed description of phantom fab- 
rication is given in references [28], [29]. The results reported 
here were obtained from a specific set of phantoms illustrated 
in Fig. 3. 

A. Homogeneous Phantom 

A homogeneous cylinder 88 nmi in diameter and 140 
mm long was constructed from a 5.5% by weight gelatin 
concentration. A small amount of polystyrene microspheres 
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(a) 



(b) 




UiKBT TrazBducer Amy 
(128 elements) 

(c) 

Fig. 3. Illustration of phantoms and defomiatioo used in experiments, (a) Cross sectional view of the phantom with a single hard or soft inclusion, (b) Cioss 
sectional view of the pluntom with three hard inclusions of different size. In boA cases the inclusions extend the entire length of the phantom, (c) Longitudinal 
view of the phantom widi a short 10 mm long cylindrical hard inclusion at the center. Four dilTerent imaging planes dirough the phantoms are also shown. 



was added to act as ultrasonic scattering centers. This scatterer 
concentration remained the same for all phantom materials. 



B. Phantom with a Single Hard Inclusion 

First, a homogeneous phantom was constructed from a 5.5% 
by weight gelatin concentration. Then, a circular, longitudinal 
hole 30 nrni in diameter was made in the center of the 
phantom. This hole was backfilled with 12% gelatin concen- 
tration producing a hard cylindrical inclusion extending the 
entire length of the phantom. The Young's modulus of the 
inclusion was estimated to be about three times larger than 
that of surrounding material. This estimate was based on an 
approximately linear dependence of the Young's modulus on 
gel concentration in this range [36]. 



C. Phantom with a Single Soft Inclusion 

This phantom was constructed similar to previous phantoms 
with the only difference that 12% gelatin concentration was 
used for the surrounding material and 5.5% gelatin for a soft 
inclusion. Again, the 30 mm in diameter soft inclusion ex- 
tended the entire length of the phantom, and the Young*s mod- 
ulus of the inclusion was estimated to be about 2.5-3.0 times 
smaller than that of surrounding material. The cross-sectional 
view of the phantom with a single either soft or hard inclusion 
is shown in Fig. 3(a), where dimensions are also given. 

D. Phantom with Three Hard Inclusions 

Three longitudinal holes of different diameter (6.5 mm, 
9.5 nmi and 13 nun) were made in an otherwise homo- 
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geneous cylindrical gel constructed from 5.5% by weight 
gelatin concentration. All three holes were backfilled with 12% 
gelatin producing hard inclusions of different size, as shown 
in Fig. 3(b). Once again, these inclusions extended the full 
length of the phantom. 

E. Phantom with a Single Short, Hard Inclusion 

A homogeneous cylinder with the same overall dimensions 
as the homogeneous phantom was constructed from 5.5% by 
weight gelatin concentration. Then, a circular, longitudinal 
hole 30 nrni in diameter was made in the center of the 
phantom. The hole was backfilled with the same 5.5% gelatin 
concentration except for a small central part, backfilled with 
12% by weight gelatin concentration^ producing a haixl, 10 mm 
long, cylindrical inclusion in the center of an otherwise homo- 
geneous phantom. The central cross-section of the phantom 
with a single short inclusion closely corresponds to that shown 
in Fig. 3(a). A longitudinal view is presented in Fig. 3(c). 

Each phantom was placed in a water tank with cylindrical 
axis perpendicular to the axis of a 3.5 MHz, 128 channel, 1- 
D transducer array attached to the bottom of the tank. The 
phantom was centered so that the image plane approximated 
the central plane perpendicular to the longitudinal axis of 
the cylindrical phantom. The tank was filled with water 
to provide contact between the array and phantom. Simple 
siuface displacements were produced by a hydraulically driven 
piston contacting the phantom at the top center over the 
entire length of the cylinder, where movement of this piston 
was controlled by measuring ultrasonic pulse arrival time 
differences to the central array elements. This piston was a 14 
mm wide, rigid, rectangular block extending the entire length 
of the phantom. The bottom of the phantom contacted die 
tank such that it did not move during piston displacement 
The deformation at the central vertical plane of the phantom 
produced by this experimental system closely corresponds to 
the deformation geometry shown in Fig. 2 and approximates 
a plane strain state. 

As discussed previously, traditional speckle tracking breaks 
down for large (compared to an acoustic wavelength) displace- 
ments because of both out-of plane motion and decorrelation 
effects due to finite strain [28]. To overcome these limitations, 
a laige set of images were recorded with small relative 
displacements but significant total displacement from the be- 
ginning to end of the set The results presented below used a 
lai^ge set of images spanning the total vertical piston displace- 
ment with 200-300 /zm steps to quantiutively estimate the 
vertical displacement and the vertical longitudinal strain. These 
images were computed by properly accumulating differential 
displacement and strain estimates between two neighboring 
images of this large set spanning the total defoimational range 
[28]. 

As reported previously [29], in the area of the central 
vertical plane of the phantom {xs = 0), the deformation pattern 
produced by the current system will closely approximate a 
plane strain state. For this state, components ui and U2 of 
the displacement vector U are functions only of xi and 
X2, and U3 = 0. With these conditions, (12) may be used 



for boundary detection of the inclusions, and (13) should 
completely describe deformations in the imaging plane. 

After recording the set of ultrasound images and computing 
displacement and strain images from diese data, the hybrid 
reconstruction of the Young's modulus proceeds in two sub- 
sequent steps. First, inclusion boundaries must be identified to 
define a nontrivial boundary value problem. That is, the closed 
curve of constant Young's modulus in the imaging plane must 
be defined. In all cases, the stress continuity condition was 
analyzed everywhere throughout the imaging plane using (12) 
to highlight inclusion boundaries. For every vertical line, i.e. 
xi is constant, and any position X2 on this line, €22(x2 — h) 
and e22(x2 +^) were computed for a fixed step h. Using these 
definitions, the relative Young's modulus F was computed at 
every point according to (12), where F greater than 1.0 means 
that Young's modulus E{x2 + h) is greater tiian E{x2 - h) 
and F less than 1.0 means that Young's modulus E{x2 + h) 
is less than E{x2 - h). 

Also, to minimize the effect of noise on these images, the 
value of F at each pixel is multiplied by the ratio of the average 
strain magnitude (622) in a large region about the pixel to 
the average strain along the center vertical line of the image. 
As discussed in [28], strain variance is relatively independent 
of strain magnitude. Based on a simple propagation of error 
analysis, the variance in F estimated by (12) is 

<^f = (r' + i)(^) (W) 

where is the variance in the strain image. Dearly, in 
regions of small average strain £22* normalizing by the average 
strain reduces large scale fluctuations in estimates of F. This 
normalization does not alter any results in high strain regions 
and only serves to minimize the effects of noise in low strain 
areas. 

The second step is die actual reconstruction of the spa- 
tial distribution of the Young's modulus by solving (13) 
for unknown fe(X). To solve (13), both components of the 
displacement vector must be accurately measured throughout 
the entire image plane. However, only one component, in our 
experimental setup the vertical component U2{xi, X2), of the 
displacement vector is measured accurately using ultrasound 
speckle tracking [28]. Fortunately, the lateral displacement 
ui(^i) ^2) can be estimated from the vertical displacement 
U2{xu ^2) using the incompressibility condition ui,i -l-U2,2 = 
0 and measurements of the lateral component of the displace- 
ment vector along any vertical longitudinal line. The specific 
program developed to solve (13) assumes that the lateral 
displacement ui(xi, X2) is zero along the central vertical 
axis of the image. This is a good approximation for all 
phantoms studied based on the vertical central plane synunetiy 
of the experimental system. Only the phantom with three 
asymmetrically located hard inclusions [Fig. 3(b)] slightly 
violates this assumption. Nevertheless, for the top part of this 
phantom containing the smallest inclusion, this assumption 
will not be in a great error. 

After defining the boundary of the region of interest (ROI) 
and the closed contour of utuform Young's modulus, (13) is 
discretized with a 2 mm grid spacing over the ROI, where 
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Fig. 4. Normalized elastic 




and (b) phantom with a single hard inclusion. 



all spatial derivatives are approximated by finite differences. 
Original displacement images are spatially low pass filtered to 
match this grid spacing before computing finite differences and 
to ensure that the model of a continuous elasticity modulus is 
well approximated. Also, it was assumed that F12 vanishes for 
our defoimational system, where the phantom was immtersed 
in water during the experiment and, therefore, gravitational 
and buoyancy forces cancel. 

The linear set of equations resulting from the discretization 
of (13) is solved by iteration, where the error on each step 
is estimated by averaging the left-hand side of (13) over the 
entire grid using the current estimate of the Young* s modulus 
distribution. From step to step, the entire modulus distribution 
over the full grid is updated based on changes in the average 
error. 

IV. Results 

The hybrid reconstruction procedure was first tested on 
the homogeneous phantom and phantoms with a single long 
inclusion. In Fig. 4, T estimated using (12) for the homoge- 
neous phantom [Fig. 4(a)l and the phantom with a single hard 
inclusion at the center [Fig. 4(b)] are presented. A step size h 
of approximately 0.4 nam was used for both images displayed 
over a 100 mm by 100 mm area, where the transducer array 
is at the bottom and the piston is located at the top. A 
quantitative gray scale is used for these images so that a T 
of 1.25 is pure white, a T of 0.75 is pure black and a T of 
1.00 is mid-gray. Note that both images are displayed over 
the same dynamic range. There are large artifacts for both 
phantoms near the top and bottom of the image (i.e., near the 
piston and the constraining bottom plate). In these regions, 
high strain variations are formed by the botmdary conditions 
resulting in erroneous estimates of P. Excluding these areas, 
there are striking differences between the F images of these 
two different phantoms. Near the top of the hard inclusion in 
Fig. 4(b) the value of F is significantly less than 1.0, indicating 
a transition from a soft to hard region. Correspondingly, near 
the bottom of the inclusion F is significantly larger than 1.0, 



indicating a transition from hard to soft region. These are 
precisely the results expected for a single hard inclusion. 

The boundaries of the single hard inclusion can be clearly 
seen in Fig. 4(b), but are even more visible if a different 
display format, such as that of Fig. 5, is used. In Pig. 5 the 
parameter 11, defined as 



n=r-i 
n = i/F- 1 



forF > 1 
forF < 1 



is presented for the homogeneous phantom [Fig. 5(a) and (c)], 
phantom with single hard inclusion [Fig. 5(b)] and phantom 
with a single soft inclusion [Fig. 5(d)]. The first two [Fig. 5(a) 
and (b)l are the n images for the images in Hg. 4 computed 
for a step size h of about 0.4 mm. These images are subjected 
to a threshold of 0.06 (i.e., H < 0.06 is black and n > 0.06 is 
white). Again, there are major differences between these two 
images and. clearly, the vertical boundaries of the inclusion in 
Fig. 5(b) are delineated. Differing from the previous two. the 
other two images of Fig, 5 [Fig. 5(c) and (d)] compare 11 in 
the homogeneous phantom and the phantom with a single soft 
inclusion where a step size h of about 4 mm and threshold 
of 0.6 were used. Both top and bottom borders of the soft 
inclusion are highlighted in Fig. 5(d), but without the fidelity 
of Fig. S(b). The homogeneous phantom in Fig. 5(c) does not 
exhibit any boimdaries in the central area of the phantom for 
either 0.4 mm or 4 mm step sizes. The parameter 11, rather than 
r. is preferable for boundary detection since the same value 
is obtained at either hard to soft or soft to hard transitions. 
Note that each pair of images is displayed over exactly the 
same dynamic range. 

Following boundary detection, reconstruction of the 
Young's modulus for the same phantoms is shown in Fig. 6. 
The reconstruction is based on numerical solution of (13), 
where it was assumed that the Young's modulus is uniform 
along the borders of the 40 by 40 mm square positioned at 
the center of the homogeneous phantom [Fig. 6(a)], phantom 
with the single hard inclusion [Fig. 6(b)] and phantom with 
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(a) 




(b) 




Fig. 5. Threshold difrcrential elastic modulus (11) estimates . 
a single hard inclusion, (c) and (d) homogeneous phantom and 
display dynamic range. 



for two pain of {duatoms: (a) and (b) homo^neous phantom and phantom with 
phantom with a single soft inclusion. Images in each pair are displayed over the same 



the single soft inclusion [Fig. 6(c)]. This square was chosen 
since it includes the boundaries detected in Figs. 4 and 5. The 
images in Fig. 6 represent the 40 by 40 mm region, where a 
logarithmic gray scale over the range 0.5 < E/Eq < 2.1 is 
used. The gray scale was selected so that the relative Young's 
modulus of 1.0 is raid-gray, dark areas represent softer material 
and bright areas harder material. The phantoms were subjected 
to 4-5 mm uniform vertical surface deformation applied 
along the entire length of the phantom. Both hard and soft 
inclusions are cleariy visible. Moreover, changes in contrast 
for both inclusions are much greater than average error in the 
reconstruction of the homogeneous phantom. This is shown 
in Fig. 6(d). where reconstructed elasticity profiles along the 
central horizontal line are compared for the homogeneous 
phantom and phantoms with hard and soft inclusion. The 
results presented in Fig. 6 demonstrate that the reconstruction 
algorithm based on (13) reasonably matches the expected 
magnitude of relative Young's modulus in these phantoms. 

As discussed previously [29], if a plane strain state is not 
present, i.e., an inappropriate model is used, then reconstruc- 



tion based on (13) will be in enor. To experimentally test 
the magnitude of these errors for a simple model system, 
the phantom with a short hard inclusion in the center of 
an otherwise homogeneous gel was used. This phantom was 
identical to the one with a single hard inclusion except that 
the inclusion was only about 10 mm long rather than the fiill 
length on the phantom. Four different planes of the phantom 
shown in Fig. 3(c) were imaged with a uniform 5 mm siuface 
deformation, where the first imaging plane is at the center of 
the inclusion, the second plane is 3 mm displaced from the 
center of the inclusion but still within the inclusion (i.e., 2 
mm from the edge of the inclusion), the other plane is 6 mm 
displaced from the center (i.e., 1 mm outside the inclusion) 
and, finally, the last imaging plane is 10 mm ftom the center 
of the inclusion (i.e.. 5 mm outside the inclusion). The results 
of boundary detection and elasticity reconstruction in this 
phantom are presented in Fig. 7. 

The first two [Fig. 7(a) and (b)] are boundary detection 
images, where the parameter n was computed for a step size of 
0.4 mm and displayed with a threshold of 0.06 (i.e., 11 < 0.06 
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Fig. 6. Reconstructed images of the Young's modulus in the cenual plane of (a) homogeneous phantom, (b) phantom with a single hard inclusion, (c) 
phantom with a single soft inclusion, (d) Profiles of the reconstructed Young's modulus in the same phantoms along the central horizontal line are plotted 
on a semi-logarithmic scale. Phantoms were subjected to 4-S mm uniform surface displacement 



is black and IT > 0.06 is white). Fig. 7(a) corresponds to 
the first imaging plane taken at the center of the inclusion. 
Inclusion boundaries are clearly visible on this image. 

Similar results were obtained for the second imaging plane 
displaced from the center but still within the inclusion. 
Pig- 7(b) corresponds to the last imaging plane S mm outside 
of the inclusion. Inclusion boundaries were not detected 
either in this or in the previous imaging plane outside of 
the inclusion. 

Based on these images, the closed curve of constant Young*s 
modulus is defined along the border of a 40 by 40 mm 
rectangle positioned at the center of this phantom, and elas* 
ticity reconstruction of the square area is performed at all 
four planes. The results are shown in Figs. 7(c)-(0 over 
exactly the same display dynamic range used in Fig. 6. In 
Fig. 7(c)» the Young's modulus reconstruction in the central 
plane highlights the presence of the hard inclusion detected in 
Fig. 7(a). This is an expected result since a plane strain state 
is well approximated at this plane, and the image reasonably 
approximates that of the phantom with long hard inclusion. 
In the second closest plane [Fig. 7(d)] the inclusion is also 
detected by the reconstruction algorithm, although artifacts 
are slightiy larger than on the previous image. Nevertheless, 
as soon as the imaging plane is outside of the inclusion. 



the inclusion is not present on reconstructed images of the 
Young's modulus distribution [Figs. 7(e) and (01. Therefore, 
within 1-2 mm of the edge, images within and outside the 
inclusion [Figs. 7(d) and (e)) are recognizable although artifact 
levels are certainly higher. At only 5 mm beyond the edge 
of the inclusion, the reconstruction looks very much like the 
homogeneous phantom. Also, the level of artifacts is reduced 
for the last imaging plane [Fig. 7(f)] compared to the previous 
image [Fig. 7(e)] since the inclusion has less influence on 
displacement/strain fields. Note that the resolution in 0:3 (i.e., 
slice thickness) of the strain imaging system is about 1.5-2 mm 
at the center of the inclusion. These results demonstrate that 
if uniform surface loading is used for phantom deformation, 
then reconstmction based on a model of a plane strain state 
[i.e., (13)] leads to reasonable results even if inhomogeneities 
are not infinite. 

Finally, images of F and 11 for the phantom with three 
hard inclusions of different sizes are presented in Fig. 8(a) 
and (b), respectively. In these images precisely the same 
parameters and display formats as in Fig. 4 and Fig. 5(a) and 
(b), respectively, were used. These images were computed 
for a step size h = 0.4 mm to highlight the fidelity of 
inclusion boundaries. All three inclusions are detected, where 
the distance between vertical boundaries closely approximates 
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the actual physical dimensions of each inclusion. The recon- 
struction of the Young's modulus distribution of the smallest 
inclusion is presented in Fig. 8(c) for a 20 mm by 20 mm 
square with inclusion approximately in the center. The same 
logarithmic gray scale as in Fig. 6 (0.5< E/E^ <1A) is used. 
Again, the size of the detected inclusion closely corresponds 
to the actual size. However, even though the inclusion is 
noticeable in this image, the level of artifacts is comparable 



with the magnitude of the inclusion's Young's modulus. Note 
that the spatial resolution of strain and displacement images is 
about half the diameter of this smallest inclusion. 

V. Discussion 

Reconstructive elasticity imaging is possible if internal 
displacement and strain fields are accurately measured with 
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(c) 

Fig. S. Results of the hybrid Young's modulus reconstruction for the phantom with three hard inclusions of different size: (a) nonnalized elastic modulus 
(T) estimate for a central cross sectional plane of the phantom, (b) threshold differential elastic modulus (11) estimate, (c) reconstructed Young's modulus 
distribution for a 20 by 20 mm square positioned at the top of the phantom. This area includes a 6.5 mm diameter inhomogencity located approximately 
in the center of the image. 



subsequent reconstniction of the elastic (Young* s) modulus. 
In this report, we have shown that longitudinal displacement 
and strain images in the limit of large (compared to an 
acoustic wavelength) surface deformation can be used to detect 
inclusion boundaries and reconstruct the Young's modulus 
distribution. 

Boundaries of local inhomogeneities can be detected even 
in cases where the simple assumptions leading to (12) are 
not strictly valid. If a plane strain state can be q)proximated, 
then accurate estimates of F will be obtained from images 
of a single strain component at boundary points where the 
components of the normal vector are equal in magnitude, i.e., 
(m)^ = ("2)^- Even at boundary points not satisfying this 
condition, either T or II will still differ fh>m unity although 
the exact value will be incorrect if the magnitude of the 
longitudinal strain is significant. At positions where either 
the magnitude of the shear strain is greater than that of the 
longitudinal strain or the longitudinal strain signal-to-noise 
(SNR) is poor, (12) cannot be used for boundary detection. 
Nevertheless, as demonstrated by measurements on simple gel- 
based phantoms, bounded inclusions can be readily detected 
with (12) if the longimdinai strain SNR is large. 



The parameters P or 11 are computed with a fixed step. 
For the results presented in Figs. 5(a) and (b), a step of 0.4 
nun was used. This size is considerably less than the 3 mm 
spatial resolution of the original strain images [28], [29]. Such 
a small step size was used to minimize the apparent extent of 
inclusion boundaries. An unfortunate consequence is that both 
r and n are severely underestimated. To gain more accurate 
measures of the relative Young's modulus at the boundary, a 
step size comparable to the spatial resolution must be chosen. 
In Fig. 5(d) the II image for the phantom widi a single soft 
inclusion is illustrated for a step size of 4 mm and a threshold 
such that n < 0.6 is black and 11 > 0.6 is white. Both top and 
bottom borders of the inclusion are highlighted, but without 
the fideHty of Fig. 5(b). Nevertheless, peak values of 11 in 
the region where (ni)^ = (112)^ are about 1.3, close to the 
value of n = 1.5 (i.e., P = 2.5) expected for this inclusion. 
A comparison between Fig. 5(b) and (d) demonstrates the 
need to compute n at several different step sizes. Small steps 
can delineate the boundaries of local abnormalities whereas 
large steps can quantify the Young's modulus along regions 
of the boundary where (ni)^ « (112)^. If the inclusion itself 
is homogeneous, then elasticity reconstruction is complete 
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since F or 11 is estimated in boundary areas where (nj)^ n 

The Young* s modulus at the inclusion boundary was cal- 
culated from a single accurately measured component of the 
strain tensor, where this estimate is only valid at points where 
(m)^ = (7^2)^ • However, to fully characterize defonnations 
of an elastic body, the complete displacement vector and 
strain tensor must be computed. For a plane strain state, 
only lateral and longitudinal components of the displacement 
vector must be measured since the other component vanishes. 
Consequently, only one longitudinal and one shear component 
of the strain tensor need be measured to represent the strain 
tensor for an incompressible material under a plane strain state. 
Parameters F or II can be more accurately computed at all 
boundary points if every component of the strain (i.e., longitu- 
dinal and shear) is measured and the more general (1 1) is used. 

If the inclusion is not homogeneous, boundary detection 
can aid more complete reconstruction. Once the boundary is 
determined and the modulus on that boundary estimated, more 
traditional numerical algorithms can be used to reconstruct 
the elasticity distribution within the inclusion. However, the 
specific approach used to solve second order partial differential 
equations of the type given in (13) depends on the relative 
magnitudes of the coefficients in the equation [40]. This 
equation generally is hyperbolic, although a degenerate (i.e., 
parabolic) case is possible. For an incompressible material, 
the degenerate case for the plane strain state corresponds 
to an in-plane translation with rotation of the volume as 
a rigid body. This means that the distribution of relative 
Young^s modulus k(X) cannot be reconstructed in this case 
without knowledge of the pressure distribution. Nevertheless, 
if the pressure distribution can be measured, then elasticity 
reconstruction can be performed directiy. 

Numerical solution of (13) for a nondegenerate case requires 
only knowledge of the elastic modulus along some boundary, 
but not necessarily the boundary of the object. Usitig inclusion 
boundaries determined by (12). the Young's modulus must 
be specified for any closed curve outside this inclusion. In 
general, however, Goursat conditions are preferable, where the 
Young's modulus need be specified only along two intersecting 
characteristic curves [40]. These characteristic curves reduce 
simply to two intersecting edges of a square region of interest 
if the ratio of shear strain to longitudinal is small over the ROI. 
Consequently, if the external deformation can be conU-olled to 
minimize the ratio of shear to longitudinal strain components 
in a square ROI, then the boundary conditions arc greatly 
simplified. Future studies will explore Goursat conditions for 
geometries more closely approximating clinical imaging. 

The quality of Young's modulus reconstruction is ultimately 
limited by the spatial resolution and accuracy of measured 
displacements and strains. Previously, a method of tracking 
relatively large displacements was presented showing signifi- 
cant signal-to-noise (SNR) improvement in both displacement 
and strain estimates [28]. These SNR improvements should 
enhance contrast resolution in elasticity reconstruction. That 
is. even if linear elasticity does not strictly hold for large 
scale deformations, Young's modulus reconstruction based on 
(12) and (13) may still be more accurate. To test this, the 



entire set of experiments on the homogeneous phantom and 
phantoms with a single hard or soft inclusion was repeated 
with a uniform surface displacement of 13.2 mm, representing 
nearly a factor of 3 increase in the external deformation. 
Reconstructions based on these measurements are presented in 
Fig. 9 over the same display dynamic range used for Fig. 6. 
Again, Fig. 9(a) is the homogeneous phantom, Fig. 9(b) is 
the phantom with a single hard inclusion and Fig. 9(c) is 
the phantom with a single soft inclusion. A comparison 
of the results in Figs. 6 and 9 clearly demonstrates that 
contrast resolution in elasticity reconstructions is improved if 
larger external deformations are used. For more quantitative 
comparison, the profiles of the Young's modulus distribution 
along the central horizontal line of Figs. 6(a) and 9(a) (i.e., 
homogeneous phantom) are presented in Fig. 9(d). Significant 
reduction in Young's modulus variations due to enhanced SNR 
is noticeable in this graph. Therefore, SNR improvement in 
strain and displacement measurements directly translates into 
enhanced contrast resolution in elasticity reconstructions. To 
optimize both spatial and contrast resolution of elasticity im- 
ages, deformational imaging methods permitting larger internal 
strains must be used for all applications. 

In the present study, only one longitudinal component of the 
displacement vector was accurately measured in the imaging 
plane. However, the lateral ui displacement can be estimated 
directly from the measiu-ed longitudinal U2 displacement using 
the principle of incompressibility, i.e., ui,i+W2,2= 0- With 
this condition, the lateral component was estimated assuming 
that it was zero along the central vertical line of the phantom. 
This is a good approximation for the experimental system 
and almost all phantoms used in the present study. In more 
general cases, however, this approximation may not hold and, 
therefore, accurate measurements of the lateral component of 
the displacement vector at most along any single vertical 
line are needed. In general, access to the lateral component 
of the displacement vector will significantly improve both 
reconstruction and boundary detection algorithms. 

Even if strain images are intimately linked to the Young's 
modulus distribution, artifacts are present in these inuiges 
due to the details of global boundary conditions, such as the 
geometry of the object, shape changes during deformation, 
specific types of external loads, etc. Reconstructive elastic- 
ity imaging has the potential to remove these artifacts and, 
therefore, produce images relatively independent of global 
boundary conditions. This is reasonable since no knowledge 
of the global boundary conditions is required to perform 
elasticity reconstruction in a selected ROI since linear elas- 
ticity specifies that true mechanical properties of a medium 
are not load/geometry dependent. At the same time, it has 
been demonstrated that the SNR in strain and displacement 
measurements has clear impact on elasticity images. Also, 
manipulations with global boundary conditions will permit 
proper control of the deformed state and. therefore, optimal 
balance between longitudinal and shear strain in the ROI. That 
is, careful control and optimization of the boundary conditions 
are necessary for accurate Young's modulus reconstruction 
even though reconstructions are theoretically independent of 
global boundary conditions. 
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Fig. 9. Reconstructed images of the Young's modulus in the central plane of (a) homogeneous phantom, (b) phantom with a single hard inclusion, (c) 
phantom with a single soft inclusion. Phantoms were subjected to 13.2 mm uniform surface displacement, representing nearly a factor of 3 increase over the 
deformation used to produce Fig. 6. (d) Comparison of the reconstructed Young's modulus magnitude profiles in the homogeneous phantoms subjected to 4 
and 13.2 mm surface displacements. These profiles along the central horizontal line of the conesponding image are plotted on a semi-logarithmic scale. 



Finally, the results presented here suggest that quantitative 
reconstruction of the elastic modulus may be possible for 
complex objects such as the human body. Full reconstruction 
of the elastic properties of soft tissue without any assumptions 
should be based on reasonably good measuFements of the 
complete 3-D spatial distribution of all necessary components 
of the displacement and strain. In the meantime, if limited 
measurements are available, then both a correct model of 
inhomogeneities and proper control of external deformation 
must be used for accurate reconstruction. 

Appendix A 
Compressible Medium 

In a compressible medium, the Poisson's ratio i/ does not 
equal 0.5 and, therefore, the components of the stress tensor 
in an isotropic, continuous, compressible medium under static 
deformation are: 

aij = ke5ij-\'2fieij (A-1) 

where 0 = divV = en + £22 + £33 is the trace of the strain 
tensor, the parameters A and /i are Lame coefficients, and e,^ 
is one component of the strain tensor defined in (4). Note die 
differences between the stress-strain relations in (2) and (A- 



1). Assuming that the Poisson's ratio depends only slightly on 
position, and defining the following parameters 

t/ 1 

'=(i+i/)(i-2i.)' "*=2(rTiJ) ^^'^^ 

the Lame coefficients can be written in the form 

A = /£;oik(X), ti = mEQk{X) 

where Eq is constant and k{X) is arbitrary and not generally 
a continuous function of position. 

First we will consider a bounded inclusion with boundary 
G. Assigning indexes **mf * and "cxt" to internal and exter- 
nal variables, respectively, at boundary points, and substituting 
r for the ratio of Young's moduli ^^^s 
continuity condition (7) becomes [26]: 

ri/e^'*n2 + 2m(£ff + eg^nz -\- eg'na)] 
= /e'"*n2 + 2m(e\?'ui + 42*^2 + ^iV^a) 

r[/e*«*n3 + 2m{e\l*ni + eg*n2 + cg'na)] 

= /e'"*n3 + 2m(€^ni + eg^na + cS^na). (A-3) 
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Note the similarity of (A-3) with (8) for an incompressible 
material. Again, if all components of the strain tensor are 
known, then any equation in (A-3) can be used to determine F 
at the inclusion boundary, and the other two equations in (A-3) 
can be used to estimate the accuracy of strain measurements 
in the neighborhood of the boundary point. However, if the 
value of Poisson*s ratio is unknown, i.e., parameters / and 
m cannot be computed and must be eliminated from the 
system of equations (A-3), the resulting equations will be 
the same as in (9). This is an expected result since in both 
compressible and incompressible media the same A6 product 
must be eliminated. Nevertheless, if the elastic properties do 
not exhibit lai^ge fluctuations inside the inclusion, elasticity 
reconstruction is complete since the value of k{X) is found 
at the inclusion boundary. 

If the function k(X.) is continuous with well defined spatial 
derivatives, then complete reconstruction of the elasticity 
within the region of interest must be performed. For this case, 
the equilibrium condition (1) can be rewritten in the form [26] 

Eok{lGii +2meii,i +2mei2,2 +2m€i3,3 ) 

+ Eoiiekyi +2mffiiA;,i +2niei2fc,2 +2mei3A;,3 ) + /i = 0 
Eok{lB,2 +2m£i2a +2me22,2 -f 2m€:23,3 ) 
+ EQ(iek,2 +2m£:i2fc,i +2m£22A;,2 -f 2me23A:,3 ) + /2 = 0) 
Eok{lS,3 +2mci3,i +2»7ie23i2 +2m£33,3 ) 

+ Eo{lSk,z +2mei3fe,i -l-2mc23*;,2 +2171633^,3 ) 
+ /3 = 0. (A-4) 

Here fte lower index after a comma means (lifferentiation with 
respect to the corresponding spatial coordinate. 

It is important to note that each equation in (A-4) is 
only a first order partial differential equation with respect to 
k{X) and, therefore, generally more preferable. However, the 
specific equations presented here for a compressible medium 
are very difficult to solve if in fact the medium is nearly 
incompressible. In this case, 6 = divXJ is small, but the 
parameter / is large (since Poisson's ratio u is very close to 0.5 
for a slightly compressible medium). Consequently, the effect 
of very small errors in estimating 8 from deformation mea- 
surements can be greatly magnified by multiplication with /. 

Similar to an incompressible material, parameter /8can be 
eliminated from (A-4), and the system of equations (10) will 
be obtained with the following changes in notation: 

^i2=(/2,i-/i,2)/m^;o, 

F22 = (/a ,2 -/2,3 )/mEQ. 

Each equation in (10), however, is a second order partial differ- 
ential equation with respect to A:(X), and requires knowledge 
of higher strain and displacement derivatives than in (A-4). 
Since these components are measured with finite signal-to- 
noise ratio, error due to higher derivatives may increase. 

If all components of the displacement and strain are avail- 
able, hybrid reconstruction of the Young's modulus proceeds 



as outlined before. First, strain and displacement images are 
processed to highlight boundaries between regions of different 
elastic modulus. This is based on either (A-3) or (9). If the 
inclusion is clearly bounded, then elasticity reconstruction is 
complete. Spatial variations of the elasticity inside the region 
of interest (ROI) are computed following boundary detection. 
After low pass filtering both displacement and strain data 
within the ROI to ensure a continuous elasticity distribution, 
die reconstruction is performed by numerically solving either 
(A-4) or (10). This will conclude elasticity reconstruction 
within that ROI. 

Again, a specific type of applied defonnation can be consid- 
ered requiring displacement and strain data only in the imaging 
(xij 0) plane. For a plane strain state, the last equation in 
(A-3) is satisfied since 713 = 0 for this case, and the first two 
equations in (A-3) reduce to 

/e«*ni -f 2m(£ff«ni eff nj) 
^ IQint^^ + 2m(e\fn, + 42^712) 
/e«*n2 + 2m(e5f m + elfui) ' 

Also, if the parameter must be eliminated from (A-5), then 
(11) will be obtained. Both (A-5) and (11) should be satisfied 
along the boundary of a clearly bounded inclusion. Note that 
these equations can also be used for boundary detection. 

For a plane strain state, (A-4) reduces to two nontrivial 
equations: 

EQk{lQ,i -|-2m£ii,i +2m€i2,2 ) 

+ Eo{iek,i -\-2menk,i +2mei2A;,2 ) + /i = 0 

Eok{lB,2 +2m£i2,l +277l£22,2 ) 

' + EQ{iek,2 +2m£ 12^,1 +2me22A:,2 ) + /2 = 0) 

/3 = 0. »(A-6) 

For the case of a slightly compressible media, (10), with 
notations for compressible media, reduces to one nontrivial 
equation similar to (13). Again, two-dimensional in-plane 
displacement/strain measurements are needed to reconstruct 
the Young's modulus in the plane X3 = 0 if either (A-6) or 
(13) is used. 

In conclusion, incompressible and slightly compressible 
media must be treated in a similar way, since strain and 
displacement components are measured with finite accuracy 
and, consequendy, small error in these measurements may 
have large impact on the elasticity reconstruction. However, 
if the media is compressible, then equations with lower order 
derivatives can be used. 



Appendix B 
simplihcatlons from spatial symmetries 

Along widi simplifications of a general model for elasticity 
reconstruction due to careful control of the deformed state, 
other simplifications from spatial symmetries are possible. 
Below, we present several types of 3-D inclusions possessing 
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simple symmetries where the problem reduces to a condition 
similar to (11) in the particular plane. For these symmetric 
inclusions, however, many points satisfy the condition na = 0, 
similar to a plane strain state. In all cases considered here, it 
was assumed that strain and displacement data are measured 
in the imaging plane (xi, X2y 0)» i.e. = 0. 

If the inclusion has a plane of symmetry, i.e.. there is a plane 
(xi, X2i 0) crossing the inclusion such that along the boundary 
the condition na = 0 is fulfilled, then at every point on that 
boundary, with no assumption about the deformed state of the 
volume V, the first equation of (9) leads to (1 1), repeated here 
for convenience: 

nin2(eff - elf) + [(nj)^ - (ui)2]eff ' 

In this equation, ni and n2 are two orthogonal components of 
the normal vector n = (ni, n2, 0) to the inclusion boundary 
in plane (xi, X2), and all terms in this equation can be directly 
obtained from 2-D measurements in the imaging plane. If 
the inclusion changes sh^ slightly along the X3 direction 
making the ns component nonzero but small in some regions 
of the imaging plane, then terms in the first equation of (9) 
containing 713 are also small such that (11) still accurately 
models the stress continuity requirements for these regions. 
This is, for example, why boundary detection based on (12) 
for the phantom with a short hard inclusion highlights the 
borders of such an inclusion even if the plane strain state 
approximation is weak. The ultimate accuracy depends on the 
relative magnimde of the out-of-plane shear strain £13 and 623 
with respect to in-plane en, and £12 strain components. 

For an inclusion with planar symmetry, the spatial dis- 
tribution of Young's modulus should satisfy the condition 
k(xi, X2, xa) = k{xi, xj, -xy}]. Then, the first equation of 
(10) for the plane (xu ^2i 0) takes the form: 

2ei2(Aj,ll -*,22 ) + 2(U2i2 )kyi2 
+ (A3U2 + tA2,33 +Wi2,l )k,i -(Asttl + t4i,33 
- t*'l2,2 )kf2 

+ A3u;i2Jfc + wi2,33 A; + F12 = 0 (B-l) 

where 

A3=aV(ax?)+aV(ax^), 

Fi2=(3/£o)(/2,l -/l,2). 

The out-of-plane displacement component U3 has been 
eliminated in this equation, and other components of the 
displacement/strain can be computed from 2-D measurements 
in the imaging plane. However, terms 1*1,33 , U2,33 and C4;i2,33 
cannot be computed from measurements in a single imaging 
plane only, and create reconstruction error if ignored. Nev- 
ertheless, if measurements are possible in three close parallel 
imaging planes, then second order derivatives can be found for 
a central imaging plane using a finite difference approximation 
and, consequently, included in the reconstruction algorithm. 



In general, however, the magnitude of these terms can be 
minimized by carefully controlling the deformed state of the 
medium. In particular, the plane strain state is one such 
controlled deformation. 

If the inclusion has axial symmetry, for example an xi-axis 
symmetrical inhomogeneity, where relative Young's modulus 
Ar(X) and inclusion boundary G take the forms 

A:(xi, X2» X3) = A:(xi, r), 
G:p(xi,r) = 0, 
r^ = {x2f'^{xz)\ 

then the first equation of (9) for the imaging plane X3 = 0 
will lead to 

» _ rxagg g,r (eir " 42^) + [(^2j?,r - (rg,i )^]e\f 
rx29a 9.r {e\f - eg*) + [(x2P,r f - (rpa ^elV ' 

Again, all terms in this equation can be directly computed 
from measured strain components in the imaging plane. Note 
that no assumption is made about the deformed state. 

For actual reconstruction of the Young's modulus k{X), 
the first equation of (10) in the imaging plane X3 = 0 and, 
therefore, r = X2, leads to an expression identical to (B- 
1) permitting measurements in three close parallel planes to 
access out-of-plane derivatives. Similar results can be obtained 
for an inclusion with X2-axis of symmetry, i.e., 

k{xu X2, X3) = A:(x2, r), 
G:9{X2, r) = 0, 

Finally, a spherically symmetric inhomogeneity is a partic- 
ular case of the axially symmetric inclusion considered above 
with spatial distribution of Young's modulus 

A:(xi, X2, X3) = k(r), = (xi)^ -f- {x2)^ + (^3)^ 

Here, for an arbitrary deformed state of volume V, the first 
equation of (9) for the plane X3 = 0 reduces to a second order 
linear ordinary differential equation rather than a second order 
partial differential equation (B-l) 

^ ^k.rr-^k.r ^ [(x] - xl)£i2 + XiX2(U2.2 -«1,1 )1 
+ - k,r (a:i(A3U2 + U2,33 +t^l2,l ) 

r 

- X2(A3ttl + t«i,33 -«*^12,2 )1 

+ Aj{ A3W12 + u;i2 ,33 ) + i^i2 = 0 (B-2) 
which has a simple form for xi = r, X2 = 0: 

2€i2 (^k,rr -^k,r^ + fc,r {A3U2 + U2,33 +^12,1 ) 
+ fc(A3U;i2 + t«;i2,33 ) + ^12 = 0. 

Again, all terms containing second order partial derivatives 
with respect to X3 can either be obtained from measurements of 
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corresponding displacement/strain components in three close 
parallel planes, or ignored if the deformed state is controlled 
SO that these terms vanish. 

In practice, however, neither a perfect plane strain state nor 
symmetric inhoraogeneities may exist. FuU reconstruction of 
the elastic properties of soft tissue without any assumptions 
should be based on complete 3-D measurements of all nec- 
essary displacement and strain components. If only limited 
measurements are available, then both a correct model of 
inhomogeneities and careful control of external deformations 
must be used to achieve a reasonably accurate reconstruc- 
tion. 
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